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We study historical dynamics of joint equilibrium distribution of stock returns in the U.S. stock 
market using the Boltzmann distribution model being parametrized by external helds and pairwise 
couplings. Within Boltzmann learning framework for statistical inference, we analyze historical 
behavior of the parameters inferred using exact and approximate learning algorithms. Since the 
model and inference methods require use of binary variables, effect of this mapping of continuous 
returns to the discrete domain is studied. The presented analysis shows that binarization preserves 
market correlation structure. Properties of distributions of external fields and couplings as well as 
industry sector clustering structure are studied for different historical dates and moving window 
sizes. We found that a heavy positive tail in the distribution of couplings is responsible for the 
sparse market clustering structure. We also show that discrepancies between the model parameters 
might be used as a precursor of financial instabilities. 
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I. INTRODUCTION 

Price formation on a financial market is a complex 
problem: It reflects opinion of investors about true value 
of the asset in question, policies of the producers, exter¬ 
nal regulation and many other factors. Given the big 
number of factors influencing price, many of which un¬ 
known to us, describing price formation essentially re¬ 
quires probabilistic approaches. In the last decades, syn¬ 
ergy of methods from various scientific areas has opened 
new horizons in understanding the mechanisms that un¬ 
derlie related problems. One of the popular approaches is 
to consider a financial market as a complex system, where 
not only a great number of constituents plays crucial role 
but also non-trivial interaction properties between them 
mm- For example, related interdisciplinary studies of 
complex financial systems have revealed their enhanced 
sensitivity to fluctuations and external factors near crit¬ 
ical events HU with overall change of internal structure 
OH]. This can be complemented by the research devoted 
to equilibrium [SHE] and non-equilibrium [laiii] phase 
transitions. 

In general, statistical modeling of the state space of 
a complex system requires writing down the probability 
distribution over this space using real data. In a sim¬ 
ple version of modeling, the probability of an observable 
configuration (state of a system) described by a vector of 
variables s can be given in the exponential form 

p(s) =Z"iexp{-^H(s)}, (1) 
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where 1-L is the Hamiltonian of a system, [3 is inverse tern- 
perature (further /? = 1 is assumed) and Z is a statis¬ 
tical sum. Physical meaning of the model’s components 
depends on the context and, for instance, in the case of fi¬ 
nancial systems, s can represent a vector of stock returns 
and TL can be interpreted as the inverse utility function 
m- Generally, TL has parameters defined by its series ex¬ 
pansion in s. Basing on the maximum entropy principle 
[min], expansion up to the quadratic terms is usually 
used, leading to the pairwise interaction models. In the 
equilibrium case, the Hamiltonian has form 

'H{s) = — h'''s — s'''Js, (2) 

where h is a vector of size N of external fields and J 
is a symmetric N x N matrix of couplings (j denotes 
transpose). The model may also involve hidden states 
(nodes), which are not directly observable, but here we 
restrict ourselves to considering the visible nodes case 
only. 

The energy-based models represented by Eq. 0 play 
essential role not only in statistical physics but also in 
neuroscience (models of neural networks [TSUin]) and ma¬ 
chine learning (also known as Boltzmann machines |20)). 
Recently, applications of the pairwise interaction models 
to financial markets have been also explored miEIHls]. 
Given topological similarities between neural and finan¬ 
cial networks (53], these systems can be considered as 
examples of complex adaptive systems |25j . which are 
characterized by the adaptation ability to changing envi¬ 
ronment, trying to stay in equilibrium with it. From this 
point of view, market structural properties, e.g. cluster¬ 
ing and networks [26ll28j . play important role for mod¬ 
eling of stock prices. Adaptation (or learning) in these 
systems implies change of the parameters of PL as finan- 
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dal and economic systems evolve. Using statistical in¬ 
ference for the model’s parameters, the main goal is to 
have a model capable of reproducing statistical observ¬ 
ables based on time series for a particular historical pe¬ 
riod. In the pairwise case, the objective is 


('^2)data ('^i)niodel; /o\ 

data — models 

where j = 1,..., iV and angular brackets denote statisti¬ 
cal averaging over time. 

Having specified general mathematical model, one can 
also discuss similarities between financial and infinite- 
range magnetic systems in terms of phenomena related, 
e.g. extensivity, order parameters and phase transitions, 
etc. These features can be captured even in the simpli¬ 
fied case, when Si is a binary variable taking only two 
discrete values. Effect of the mapping to a binarized sys¬ 
tem, when the values Si = -|-1 and Si = —1 correspond to 
profit and loss respectively, is also studied in the paper. 
In this case, diagonal elements of the coupling matrix, 
Jii, are zero because sf = I terms do not contribute 
to the Hamiltonian. It is worth stressing that the cur¬ 
rent investigation develops ideas outlined in the previous 
studies [l5l[2Tli2^ in a way that the effect of binarization, 
comparison of learning algorithms, evolution and scaling 
properties of the parameters distributions are studied in 
a more systematic fashion. 

The paper is organized as follows. In Section]^ basic 
statistical definitions and inference methods for h and 
J are presented. In Section m effect of binarization 
of stock returns and historical evolution of the model 
parameters are discussed. Finally, the main findings of 
our investigation are summarized in Section [IV| 


II. DATA AND METHODS 



1996 2000 2004 2008 2012 

(b) DFT of s (Amplitude) 


0.04 - 
0.02 




0.00 

0.04 

0.02 


II i\ .liilLlliJllijuiiLiliiLlhilLiiyili 




0.80 


0.00 






0.0 


0.1 


0.2 


0.3 


0.4 


0.5 


(jj (days'^) 


FIG. 1. Historical dynamics of the mean raw and binary 
returns (a) and amplitudes of their Discrete Fourier trans¬ 
forms (b). The distance between two labeled dates is 1000 
trading days. A few major financial crises are highlighted 
with the light green background: (1) Asian and Russian cri¬ 
sis of 1997-1998, (2) dot-com bubble, (3) U.S. stock market 
downturn of 2002, (4) U.S. housing bubble, (5) bankruptcy of 
Lehman Brothers followed by the global financial crisis, (6) 
European sovereign debt crisis. The number in each panel 
shows overall historical correlation between the correspond¬ 
ing series. Historical values of average return, economic cycles 
and frequency of crashes are preserved for the binary returns 
despite the maximum magnitude being bounded. 


equally weighted previous values of daily log-returns (in¬ 
cluding the current one). For each set of chunks, one 
can calculate different statistical characteristics, such as 
average 


We study historical dynamics of the U.S. stock mar¬ 
ket using N = 71 stock prices time series [29] from 
the Standard & Poor’s 500 index (hereafter S&P 500) 
listed in Table [l| We analyze discrete daily closing prices, 
S'i(t), starting from 1990 till 2013 (5828 trading days), 
which are converted to logarithmic returns, (t) = 


A. Basic statistical analysis of financial time series 

The top panel in Fig. [^a) shows historical data for 
the average stock return. In order to extract long-term 
trends from the time series, we employ a simple moving 
window (or simple moving average, SMA) approach. It 
acts like a low-pass filter, averaging out high frequency 
components, which are usually related to noise. Within 
the SMA approach, data is divided into chunks (or win¬ 
dows) of size T, assuming the time series to be stationary 
on this scale. In this case, tth chunk corresponds to T 


= (4) 

t'=0 

and covariance matrix, C {N x N), with the elements 

Cij = {SiSj) — {Si){Sj), (5) 

where T > N is assumed for the covariance matrix to be 
positive definite. Hereinafter, angular brackets ( ) de¬ 

note averaging over time (historical values), while bar 
denotes averaging over index (vector or matrix elements). 
It is also possible to investigate nonlinear dependence 
between time series using more sophisticated statistical 
concepts [30] or nonlinear data transformations, however 
only the simplest linear case is considered in the current 
paper. Series variance is autocovariance, erf = Ca, where 
a denotes standard deviation or volatility in finance. Usu¬ 
ally, SMA volatility serves as the simplest risk measure 
quantifying stability of returns. In order to quantify de¬ 
viation from the normal distribution, it is also useful to 
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define higher-order moments, such as skewness 


B. Equilibrium Boltzmann learning methods 


(Skew (si)) 



- (s*) 


O'* 


and kurtosis (also known as excess kurtosis) 


(Kurt {si)) 



We harness inference methods based on maximizing of 
the model’s likelihood £(h, J | In the equilibrium 

case, exact learning of the Hamiltonian [Eq. ([^] parame¬ 
ters implies solving Eq. @ in a self-consistent way, where 
corrections Shi and 5Jij on each learning step can be cal¬ 
culated as 


( 7 ) 


— TJfi ((s*)data ('5*)model) j 
5Jij — rjj ((s*Sj)data (■5*5j )niodel) ■ 


( 11 ) 


which both equal zero in the Gaussian case. Indeed, SMA 
filter allows one to extract long-term trends in the market 
and, for a large value of N, various moments of the distri¬ 
bution of returns can be used to identify market crashes 
(Figs. i and[^. Henceforth, we will call the considered 
portfolio “market” as it represents a big number of the 
top companies from S&P500 index. Also, we provide cor¬ 
responding confidence intervals for these moments using 
basic bootstrapping algorithm (sampling with replace¬ 
ment) which are depicted in Figs. and However, 
a more comprehensive discussion on the moments esti¬ 
mates for the non-stationary time series, besides moving 
window and correlated variable effects [3l], would also 
require considering nontrivial long-term memory effects 
and intraday correlations |32L 133) inherent to financial 
time series, what goes beyond the scope of the current 
investigation. 

As mentioned in the Introduction, to reduce the 
amount of data required for our inference, we focus on a 
binarized version of the returns and not on the raw data. 
We thus define the binarized version of the returns as 

sf^ = sign{sD. ( 8 ) 


This procedure also mitigates the scaling problem of mul¬ 
tiple series since it assigns the same value of stock return 
independently of its absolute value. Another common 
technique to deal with this problem is standardization 


std 

■’i 


,raw /„raw\ 

•i ~ w* ) 


( 9 ) 


In this case, correlation matrix, Q {N x N), is a normal¬ 
ized covariance matrix with the elements 


Qij 


(JiCTj 


( 10 ) 


Standardization procedure is known to preserve market 
behavior and widely used for the statistical analysis of the 
financial time series m Effects of these transformations 
defined by Eqs. ®-(10l are shown in Figs. [Iip and will 


be discussed further in Section HI where some simple 
statistical properties of the binarized return versus the 
raw and standardized returns are compared. Before this, 
however, we briefly describe the inference procedure. 


Here, rjfi and rjj are learning rates, (•)data are empirical 
(observed) moments and (-(model are the moments sam¬ 
pled from the model using Monte Carlo (MC) methods. 
The exact learning algorithm always yields optimal val¬ 
ues for h and J if there are no hidden nodes in the system 

m- 

Being in general slow, the exact learning algorithm 
might be substituted by the approximate inference meth¬ 
ods |36j which are based on expansion of the free energy 
of a system for small fluctuations around its mean value. 
The first-order (naive) approximation within the mean 
field theory (nMF) gives 

jnMF=A-l-C-l, 

^ tanh-i(s*) - E 


where Aij = (1— {si)'^)Sij and Sij is the Kronecker delta. 
Here, taking into account the diagonal element Ja (which 
is usually discarded) for the calculation of correspond¬ 
ing hi improves accuracy of the approximation, being 
known as the diagonal-weight trick [Ml- The second- 
order correction to the nMF approximation requires solv¬ 
ing Thouless-Anderson-Palmer (TAP) equations 


(c-i 

/jTAP = ;,nMF 


ij 


TAP 


-2{jJ^^f{s,){s,), 


N 


i=i 


(13) 


where Eq. (12) should be used instead for the calculation 


of the external fields if the diagonal-weight trick is used. 

Another class of approximations can be derived using 
expansion of the free energy in pairwise correlations. The 
simplest independent-pair (IP) approximation assumes 
independence of every stock pair from the rest of the 
system. In this case, couplings and external fields can be 

found as m 




hr = iln 


(l+rrii-i-mj S-C,* - )(l — m;—rrij +<?** ) 

(l—mi-\-mj—Cr)(l+mi—mj—Cr) 

N 

j 


( l+mA 

l-rrii J 


(14) 


where C*j = Cij -\- mimj. Sessak and Monasson (SM) 
derived higher-order corrections to the IP approximation 
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FIG. 2. Historical dynamics of the first four temporal moments of the distribution of mean market return (a). Top-bottom: 
temporal mean, standard deviation, skewness and kurtosis calculated using SMA window of 250 days (approximately one 
trading year) for the raw blue), standardized green) and binary (s*’'", red) returns of 71 U.S. stocks (Table |^. 95% 

confidence intervals for the moments of the distribution of s'’*" are calculated using bootstrapping algorithm and denoted with 
the gray dashed lines. The number in each panel corresponds to overall historical correlation between the time series for 
and Dependence of these overall correlations on the moving window size is shown in (b). Binary returns behave similar 
to raw and standardized returns, however information about kurtosis is completely lost. 



FIG. 3. Historical dynamics of the first four moments of the distribution of off-diagonal elements of covariance (C^^", blue) 
and correlation green) matrices of raw returns, and covariance matrix of binary returns (C*^'", red) calculated using 

SMA window of 250 days (a). Top-bottom: mean, standard deviation, skewness and kurtosis of the off-diagonal elements of the 
matrices. 95% confidence intervals for the moments of distribution of are calculated using bootstrapping algorithm and 
denoted with the gray dashed lines. The number in each panel corresponds to overall historical correlation between the time 
series for and C*”". Dependence of these overall correlations on the moving window size is shown in (b). Binarization 

makes covariance matrix similar to the correlation matrix of raw returns. 


in a closed form using other terms in the perturbative 
correlation expansion m 

tSM _ rnMF i rpair__ 

hSM ^ Impair 

It is also worth noting that there have been developed a 
few other approximate inference schemes tailored for dif¬ 
ferent system regimes, such as a pseudo-maximum likeli¬ 
hood inference using all data |38j . minimum probability 


flow |39j or mean field approximations for low tempera¬ 
tures |40j . For example, pseudo-maximum likelihood in¬ 
ference [ 31 ], being more suitable for sparse connections, 
can be studied in a context of the correlation clustering 
structures described in Subsection IIII PI 
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III. RESULTS AND DISCUSSION 
A. Effect of binarization 

The mapping defined by Eq. ([^ obviously affects infor¬ 
mation contained in the time series. In order to estimate 
this effect, we compare time series for the average raw 
and binary returns first. As Fig.j^shows, correlation co¬ 
efficient between historical values of and s'’™ is very 
high (0.91). Moreover, correlation between amplitudes 
of their Fourier transforms is also high (0.80), suggest¬ 
ing that signatures of economic cycles and frequency of 
market crashes are preserved in the binarized time se¬ 
ries. However, maximum magnitude of binary returns is 
obviously bounded by the dehnition. 

Comparison of the filtered time series using SMA also 
suggests the idea that the information about the market 
trends is preserved in the binarized time series. With 
this aim, we compare historical evolution of the first four 
moments of the distribution of average binarized versus 
raw and standardized returns. The results presented in 
Fig-i indicate that the binary returns behave similarly 
to the raw and standardized returns, preserving dynam¬ 
ics of the first two moments and less so about the third 
moment, while information about kurtosis is lost for all 
T [Fig. [^b)]. Regarding the pairwise correlations. Fig¬ 
ure [3] shows that the covariance matrix of the binarized 
returns becomes similar to the correlation matrix of the 
raw returns. Indeed, their off-diagonal elements follow 
similar distributions with very high correlation between 
means. For higher-order moments, correlation decreases 
with T [Fig.j^b)]. 

Another way to explore effects of the binarization pro¬ 
cedure is to study eigenvalue distribution of correlation 
matrices for the original and transformed time series, 
since it is known to be different from the distribution of a 
random correlation matrix m- With this aim, we com¬ 
pare historical dynamic of the four biggest eigenvalues, 
which do not match the Wigner law. As Fig.j^shows, all 
four values are well preserved and the biggest one corre¬ 
sponding to the “market mode” is in remarkable agree¬ 
ment after binarization independently of moving window 
size. 

Having validated that binary returns indeed capture 
market historical behavior similarly to standardized re¬ 
turns, we proceed to the Boltzmann model inference and 
evolution of its parameters. 


B. Comparison of approximate and exact learning 
methods 

Following the SMA approach with T = 250 trading 
days (approximately one trading year), we study histori¬ 
cal evolution of the external fields and couplings inferred 
for each chunk of the binarized time series. We em¬ 
ploy four approximate (nMF, TAP, IP and SM) learning 
methods described in the previous section and compare 


them with the exact learning, where MC sampling is used 
for the latter. 

Figure]^ shows comparison of the inferred parameters 
for four different historical dates. Without use of the 
diagonal-weight trick, inference of external fields in the 
nMF case works better when the dates far from mar¬ 
ket crashes are considered (top row in the two leftmost 
columns, corresponding to 22 Dec 1997 and 09 Jan 2002). 
TAP correction slightly improves the inference accuracy 
however producing overestimated values in comparison 
with the exact learning algorithm (second row). Both 
MF approximations perform worse for the data near mar¬ 
ket crashes (two top rows in the two rightmost columns, 
corresponding to 09 Oct 2008 and 13 March 2012), while 
the diagonal-weight trick allows to achieve almost per¬ 
fect accuracy in both cases (red triangles in the first two 
rows), justifying the need of high-order corrections. The 
IP algorithm (third row) shows much worse performance 
for the all dates except 09 Jan 2002. Although higher- 
order SM corrections (third row) considerably improve 
accuracy of the IP external fields, it is still lower than in 
the MF cases. Couplings inferred using both MF approx¬ 
imations show the same trend, being more/less accurate 
far from/near crashes (fourth row). Although the bulk 
of the MF couplings is in an excellent agreement with 
the exact couplings even near crashes, positive outliers 
are overestimated. The couplings estimated using the IP 
algorithm are in a very poor agreement with the exact 
couplings for the all dates considered (bottom row), while 
the use of the SM correction yields considerable improve¬ 
ment (bottom row). Moreover, it correctly captures the 
biggest outliers, outperforming the MF algorithms. 

Historical dynamics of the inference quality for approx¬ 
imate methods is depicted in Fig. where normalized 
root mean square error 


NRMSE(x,y) 


{xj - Vi) 


(16) 


and correlation between the model parameters inferred 
using different methods are presented. It confirms that 
the diagonal-weight trick considerably improves inference 
quality of the MF external helds, while couplings are sys¬ 
tematically overestimated. Within the small correlation 
approximation, IP algorithm does not perform well for 
all historical dates, while quality of the SM couplings is 
comparable to the MF cases (however being still lower 
in general) and decreases for the periods where higher 
correlations are observed. Thus, the observed behavior 
indeed justifies use of TAP as a reliable approximation 
of true couplings and external fields (making use if the 
diagonal trick), which has been extensively used in the 
previous studies [HJIIIHISI- However, special care should 
be taken about overestimated positive outliers, where SM 
algorithm can be helpful. 

To validate results obtained by the exact learning al¬ 
gorithm, we sampled stock returns from the model using 
MC sampling. The single and pairwise empirical mo- 
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FIG. 4. Historical dynamics of the four biggest eigenvalues of the correlation matrix of raw returns (green) and covariance 
matrix of binary returns (red) calculated using SMA window of 250 days (a). The number in each panel corresponds to overall 
historical correlation between the time series. Dependence of these overall correlations on the moving window size is shown in 
(b). Binarization preserves biggest eigenvalues of the correlation matrix. 


merits are nicely recovered from the exact model (fig¬ 
ures are not shown), while the third-order covariances, 
((si - (si)) (sj - (sj)) (sfc - (sfe))), are almost not cap¬ 
tured by the model (Fig. [^. Although the inability to 
capture third-order correlations could be due to the lack 
of data for estimating the real third-order covariances, 
this is unlikely to be the case as the degree of mismatch 
between the third-order covariances of the Ising model 
and the data is the same even when somehow bigger 
amounts of data are used for their estimation. Never¬ 
theless, in spite of the fact that individual third-order 
covariances are not captured well, historical dynamics of 
their mean can be recovered from the model (Fig. [^. 


C. Distribution of inferred parameters 


Figure shows histograms of the external fields and 
couplings inferred using the exact learning algorithm for 
three different window sizes: T — 250,1000 and 5000 
trading days. The distribution of external fields is close 
to the Gaussian and does not possess outliers indepen¬ 
dently of T. The bulk of the couplings is also distributed 
normally, however a heavy positive tail is present. For 
bigger window sizes, the Gaussian bulk component of 
the inferred couplings becomes less prominent and the 
tail starts to dominate. Although this behavior has been 
previously reported in Ref. [22] , a role of the tail has re¬ 
mained unclear. We will address this question in the next 
subsection. 


In fact, as shown in Fig. 10 all moments of the dis¬ 
tribution of Jij scale with T. On the contrary, there is 
no obvious dependence for hi, which behavior is close 
to the external fields inferred on randomly shuffled time 


series (where each element is swapped with other ran¬ 
domly picked element using uniform distribution within 
the moving window chunk). Higher value of h inferred 
on randomly shuffled time series can be explained by the 
fact that it compensates positive (ferromagnetic) contri¬ 
bution to the mean return stemming from the positive 
mean of the inferred couplings: In the shuffled case, J 
becomes zero while all (s^“) remain unaffected, there¬ 
fore this effect should be compensated by the external 
fields. 

Historical dynamics of the moments of the distribu¬ 
tions shown in Fig. [^indicates that the average external 
field is strongly correlated with the mean return (0.90), 
while higher moments, being almost stable over the his¬ 
torical period considered, seem not to convey any par¬ 
ticular information about marke^j^ehavior. Without use 
of the diagonal-weight t rick , h is completely incon¬ 
sistent with h [Fig. 11 [a)]. Although behave 

more similar to they are significantly overesti¬ 

mated. External fields inferred using the IP and SM 
algorithms show similar dynamics to and 

respectively, however being even more greatly overesti¬ 
mated (figures are not shown). At the same time, the 
diagonal-weight trick allows to achieve almost perfect ac¬ 
curacy in both MF cases [Fig. [lT|[b)]. Distribution of 
couplings has stable small positive mean corresponding 
to ferromagnetic interaction, however with almost half 
of couplings being negative, i.e. the system is likely to 
exhibit frustrated configurations. The SM algorithm in 
general performs better than TAP for couplings estima¬ 
tion except for their mean value J , which becomes 


inconsistent with J for the historical periods with 
highly correlated market [Fig. c), top panel]. 

It is also interesting to note that the standard deviation 
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FIG. 5. Comparison of external fields (three top rows) and couplings (two bottom rows) inferred using the exact and approximate 
(nMF, TAP, IP and SM) learning algorithms for four different historical dates and SMA window of 250 trading days. The 
diagonal-weight trick is essential for correct inference of external fields in the mean field (nMF and TAP) cases. Both mean field 
approximations overestimate couplings with big absolute values. Accuracy of the IP algorithm is very low for both external 
fields (third row) and couplings (bottom row). SM corrections significantly improves the IP results and outperforms the TAP 
approximation for the positive outliers in the distribution of couplings. 
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FIG. 7. Comparison of the third-order covariances calculated 
using the real data and 50000 MC samples from the exact 
Ising model for two different moving window sizes. The Ising 
model is not capable of recovering observed individual third- 
order covariances independently of moving window size. 
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FIG. 6. Comparison (normalized root mean square error and 
correlation coefficient) of the external fields and couplings in¬ 
ferred using different methods: (a) nMF versus exact (blue), 
TAP versus exact (green) and nMF versus TAP (red); (b) 
IP versus exact (blue), SM versus exact (green) and IP ver¬ 
sus SM (red). Without use of the diagonal-weight trick, MF 
approximation quality of the external fields is very low and 
significantly drops during financial crashes [(a), two top pan¬ 
els], while making use of the trick improves it considerably 
[(a), two middle panels]. SM algorithm considerably improves 
quality of IP external fields, however being still lower than 
MF with the diagonal trick [(b), two top panels]. Couplings 
inferred using both MF methods are almost the same, with 
approximation quality being lower during the periods of fi¬ 
nancial crises ](a), two bottom panels]. The same behavior is 
observed for SM couplings [(b), two bottom panels). 


of couplings far from the crashes almost linearly increases 
over the whole period considered from 0.14 in 1996 to 0.2 
in 2013 [the second panel in Fig. [lT|[c)]. This observa¬ 
tion gains more meaning when we note that the standard 
deviation of the couplings in 1996 equals approximately 
to the standard deviation of the couplings inferred on 


FIG. 8. Historical evolution of the first four moments of dis¬ 
tribution of the third-order covariances calculated using real 
(blue), randomly shuffled (gray) and 5000 MC samples (for 
each historical date) from the exact Ising model (red) data. 
The exact Ising model is only capable of capturing histori¬ 
cal dynamics of the mean value of the empirical third-order 
covariances. However, all higher moments for the real data 
behave similarly to the ones calculated on randomly shuffled 
time series. 


randomly shuffled returns, while its value is almost twice 
bigger than one for the shuffled time series in 2013. Also, 
during the biggest market crashes, standard deviation 
of couplings has jumps because interconnections on the 
market become tighter as a result of the herding behavior 
during a financial turmoil: Until system is not adapted 
to a new economic reality, prices tend to move collec¬ 
tively with overall market performance as a benchmark 
[42]. The heavy positive tail, which can be characterized 
by higher order moments, increases for some historical 
periods. A reason for it is unknown at the moment. 
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FIG. 9. Histograms of the external fields and couplings in¬ 
ferred using the exact learning algorithm for 27 Jan 2010 for 
three different moving window sizes. The red curve denotes 
the Gaussian fit. For small moving window sizes, the bulk of 
couplings is distributed normally, while a heavy positive tail 
dominates for bigger sizes of moving window. 


D. Industry-related clustering structure 

It is well known that a stock market possesses a hier¬ 
archical clustering structure which can be detected using 
correlations [33 HSj or couplings between stock prices 
m or trading volumes [23] . One of the most popular 
techniques employed to find related structures is based 
on the minimum spanning tree (MST) algorithm. For 
instance, Prim’s algorithm—a basic MST construction 
algorithm—consists of the three following steps: 

1. Initialize a tree with a single stock, chosen arbitrar¬ 
ily from the all stocks. 

2. Find the stock not yet in the tree which has the 
strongest coupling (biggest value of Cij or J^-) to 
a stock in the tree and include it to the tree. 

3. Repeat step 2 until all stocks are in the tree. 

Figure [T2| shows examples of the MST constructed using 
covariance and coupling matrices of binary stock returns, 
which are similar to the previously reported in Refs. [15j 
and [33]. The remarkable feature of these structures is a 
manifestation of industry sector clustering. Here, we de¬ 
fine an industry sector cluster as a connected subset of the 
tree where all stocks belong to the same industry sector, 
i.e. each cluster member interacts to the other members 
only through the stocks from the same sector. Further, 
we denote cluster size as Nm,k, where m = 1,..., M is a 
sector index {M = 9 is the total number of sectors listed 
in Table |T| and k = 1,.. ., Km is an index of the cluster 
{Km is a number of such clusters for the sector m). 

As mentioned above, there can be many clusters for 
each industry sector. In order to estimate overall indus¬ 
try clustering degree of the market, let us introduce a 
simple metric based on finding clusters of the maximum 


size 


1 


M 


Qmst — ^ ' rnax Afj72,/c- 

iV k 

m—1 


(17) 


Intuitively speaking, a small value of Qmst corresponds 
to the case where stocks do not group with each other 
based on which industry sector they belong to. Its mini¬ 
mum value, M/N, is defined by the situation where the 
biggest cluster for each sector has only one stock, i.e. 
Km equals to the total number of stocks in the sector 
m. The maximum value of Qmst is 1, which corresponds 
to the perfect industry clustering structure when there is 
only one cluster for each industry sector {Km = !)■ This 
clustering measure shows interesting dependence on the 
size of moving window (Fig. [T^), suggesting an increas¬ 
ing degree of sectoral connectedness of the stock market 
for bigger time windows as inferred by the Ising model. 
Also, the degree of connectedness increases with devia¬ 
tion of the distribution from the Gaussian measured by 
skewness and kurtosis. 

To further investigate the network structure and clus¬ 
tering degree of J, we perform the clustering analysis 
based on two different filtering procedures, namely, con¬ 
sidering only a subset of (i) couplings Jij and (ii) eigen- 
modes of J corresponding to different eigenvalues A^. 
With this aim, we choose thresholds J*’*' and A**' and 
construct MST only using values Jy ^ and Xi ^ 


A*^ respectively. Figure 14'a) shows that the biggest 


drop/increase of Qmst occurs only when the a small per¬ 
cent of the strongest couplings is excluded/included for 
MST construction. In a similar way, it is also sensitive to 
discarding the biggest eigenvalues [Fig.[l4|(b)]. Thus, one 
might conclude that distribution of couplings indeed can 
be viewed as a mixture of the Gaussian bulk and heavy 
positive tail which contains all information about market 
clustering structure. From this perspective, use of sparse 
regularization methods for couplings inference, for exam¬ 
ple, £i-regularization discussed in Ref. jH], would be of 
practical interest. 

Finally, it is also worth noting that intraday inter¬ 
nal structure of couplings is neither stable (quenched) 
nor completely random (annealed), somehow preserving 
a clustering structure with the diameter however being 
smaller near market crashes (figures are not shown, see 
for example Ref. [ 15 ] )■ These non-trivial features of a 
coupling matrix will be studied in more detail in the fu¬ 
ture works. 


E. Scaling of inferred parameters 

In order to study extensive properties of the system, 
we investigate scaling properties of the parameter distri¬ 
butions with number of stocks, which are usually char¬ 
acterized by the scaling exponent A^“. We estimate its 
average value for each distribution moment over scaling 
exponents for randomly selected subsets of stocks of dif¬ 
ferent size. 
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FIG. 10. Scaling of the first four moments of the distribution of exact external fields (top row) and couplings (bottom row) with 
moving window size T for two different historical dates and randomly shuffled time series for 27 Jan 2010. Higher value of the 
mean value of the external fields inferred on randomly shuffled time series compensates positive (ferromagnetic) contribution 
to the mean return of the mean value of couplings, which is zero for the shuffled time series. 


(a) h - nMF - TAP -Exact 




1996 2000 2004 2008 2012 



1996 2000 2004 2008 2012 


FIG. 11. Historical dynamics of the first four moments of distribution of the external fields inferred without (a) and with (b) use 
of the diagonal-weight trick using the nMF (blue), TAP (green) and exact (red) inference algorithms; off-diagonal couplings (c) 
inferred using the SM (blue), TAP (green) and exact (red) inference algorithms. Moments of the distributions of the couplings 
and external fields inferred using the TAP algorithm on randomly shuffled returns are denoted with the gray dashed line. Both 
mean field approximations incorrectly estimate external fields without use of the diagonal-weight trick and tend to overestimate 
couplings. The SM approximation allows to correctly infer the biggest couplings, however their mean is incorrectly captured 
for the historical periods with high correlations. During periods of crisis, increase of couplings strength is observed. 


As Fig.[^(top row) shows, the distribution of external 
fields does not possess any particular scaling law, except 
the mean, which has a close to —0.75. The other mo¬ 
ments scale similar to the corresponding moments of the 
external fields inferred on randomly shuffled time series. 
Scaling properties of distribution of couplings depends 
on the size of moving window (Fig. 15 bottom row). 
When T is small, empirical couplings show similar scal¬ 
ing features to the couplings inferred on randomly shuf¬ 
fled returns. However, scaling of the mean and standard 
deviation becomes closer to the properties of the normal 
distribution with growth of T. This dependence might 
be related to the presence of finite-size effects, when use 
of a small number of historical values is not enough to 
estimate true correlations on the market. These results 


are similar to the scaling properties obtained in Ref. [52] , 
however their SMA window dependence has been shown 
for the first time. 

There are two extreme ways in which a change in the 
distribution of couplings can arise as one increases the 
size of the observed system. One is that the structure be¬ 
tween the previous stocks entirely changes by adding new 
stocks. Alternatively, the couplings could only change 
their absolute magnitude, while they maintain their mag¬ 
nitude relative to one another. To better understand 
where in this spectrum our financial market exists, we 
performed analysis similar to Ref. [19]: We chose a ran¬ 
dom subset of 20 stocks and analyzed couplings between 
them for different total number of stocks taken into ac¬ 
count for the inference (including the original 20). Fig- 
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FIG. 12. Minimum spanning tree for the covariance matrix 
(a) and corresponding exact couplings (b) for 27 Jan 2010 cal¬ 
culated using SMA window T = 4000 trading days. Similar 
industry-related clustering structure is observed in both cases. 
The considered sectors are Healthcare (red), Consumer Goods 
(blue), Basic Materials (green). Financial (cyan). Industrial 
Goods (purple). Services (yellow), Technology (orange). Con¬ 
glomerate (magenta) and Utilities (dark blue). The graphs 
are visualized using the NetworkX Python package |45| . 



FIG. 13. Quality (industry sector clustering degree) of a min¬ 
imum spanning tree, Qmst, depending on moving window size 
for exact and TAP couplings, and covariance matrix (27 Jan 
2010). The gray dashed and dotted lines denote mean and 
99.7% confidence interval respectively for the quality of MST 
built on randomly shuffled time series. Quality of MST for 
couplings increases with deviation of their distribution from 
the Gaussian, characterized by skewness and kurtosis. 


ure 16 shows that the biggest/smallest values of Jij re¬ 
main the same with growth of N and their scaling be¬ 
comes closer to the normal distribution for bigger time 
windows. This behavior also suggests that important fea¬ 
tures of market connectivity are preserved with the num¬ 
ber of stocks. 



FIG. 14. Cutoff analysis for the coupling matrix (a) and its 
eigenvalues (b) for two SMA windows T = 1000 and T = 5000 
trading days (27 Jan 2010). Histograms of couplings and cor¬ 
responding eigenvalues (left column), and quality (industry 
sector clustering degree) of a minimum spanning tree, Qmst, 
(left and right columns) depending on the positive (black 
dashed line) and negative (blue solid line) cutoffs, J*** and 
A**', where all values above or below a cutoff, respectively, are 
discarded for MST construction. The red curve denotes the 
Gaussian 6t. Discarding both biggest couplings and eigenval¬ 
ues significantly affects quality of MST, while negative cou¬ 
plings and small eigenvalues do not contribute to the market 
clustering structure. 


F. External and internal influence 

Considering the two terms in the system’s Hamiltonian 
[Eq. ([^], it is also possible to define internal and exter¬ 
nal influences in the market. For this purpose, exter¬ 
nal fields can be interpreted as the influence of external 
factors, h®’'* = h, while couplings define internal bias, 
= (s‘'')J (in the ME sense) [22]. In this case, exter¬ 
nal contribution corresponds to the individual stocks bi¬ 
ases which come from outside the market, while internal 


one is solely defined in terms of internal market interac¬ 
tions. Similarly, one can also define two energy terms as 
n = -f where (s). Fig¬ 

ure shows that both energies have almost the same 
order of magnitude over the historical period considered, 
while near the major crashes E®”^* is more than 10 times 
bigger than . The ratio between the mean biases also 
possesses interesting historical dynamics. Being in prin¬ 
ciple strongly correlated with the mean return (0.9 for 



















































































































12 



FIG. 15. Average scaling exponent with number of stocks, a, for the exact external fields (top row) and couplings (bottom 
row) distributions depending on moving window size, T, for two different historical dates and randomly shuffled time series for 
27 Jan 2010. Scaling exponent of the mean of external fields is close to —0.75 for big values of T, while the higher moments 
behave similar to the external fields inferred on randomly shuffled time series. Empirical couplings scale similar to the couplings 
inferred on randomly shuffled time series for small window sizes, while their scaling properties become closer to the properties 
of the Gaussian distribution for bigger values of T. 






0 200 400 

index 


FIG. 16. Scaling of the 380 couplings between a randomly selected subset of 20 stocks depending on total number of stocks used 
for the inference (including original 20) for 27 Jan 2010 and T = 1000 (top row) and T — 5000 (bottom row): 10 biggest/smallest 
couplings (first column), mean and standard deviation of the couplings (second column) and their visualization (third column). 
Relative magnitude of the couplings decreases as the number or stocks grow, while the underlying network structure is preserved. 


and 0.99 for ft,'”*) discrepancies between them might 
be used as a leading indicator of financial instabilities. 
Away from the periods of crisis, their ratio is almost sta¬ 
ble, while divergent behavior is observed before the U.S. 
market crashes (two bottom panels in Fig. 17). Possi¬ 
ble explanation of the observed behavior from a financial 
point of view is still an open question. 


IV. CONCLUSION 


We have investigated various aspects of application 
of the pairwise interaction model to financial time se¬ 
ries. The model, being parametrized by external fields 
and couplings, is used for the approximation of the joint 
equilibrium distribution of stock returns. Since the con- 
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1996 2000 2004 2008 2012 


FIG. 17. Comparison of external and internal biases for the exact Ising model. Top-bottom: ratio of external and internal 
energies contributing to the Hamiltonian, absolute value of the ratio of the mean external and internal biases and its sign. 
Major market crashes are preceded by growth of the external energy and discrepancy between the biases. 


sidered learning algorithms require use of binary vari¬ 
ables, the logarithmic returns are binarized using the 
sign function. Effect of the binarization suggests that the 
distribution of binary returns captures the distributions 
of raw and standardized returns to a good degree, pre¬ 
serving information about market correlation structure, 
economic cycles and market crashes as well as industry- 
related market clustering structure. 

The model parameters are inferred using approximate 
and exact learning algorithms. Mean field approxima¬ 
tions nicely recover bulk of the couplings except outliers, 
which can be correctly inferred using small correlation ex¬ 
pansions. External fields are in the almost perfect agree¬ 
ment with the exact algorithm if the diagonal-weight 
trick is used. The obtained results also suggest that the 
quality of approximate inference methods drops in the 
periods of financial crises due to increase of magnitude of 
the correlations observed. For the historical period con¬ 
sidered, the distribution of external fields is close to the 
Gaussian and does not possess any outliers independently 
of moving window size. On the contrary, the distribution 
of couplings possesses a heavy positive tail, which starts 
to dominate over the Gaussian bulk for bigger moving 
window sizes. Mean of external fields decreases with the 
number of stocks while their standard deviation remains 
almost constant, corresponding to the standard deviation 
of the external fields inferred on randomly shuffled time 
series. Scaling properties of the distribution of couplings 
depends on the moving window size, becoming closer 
to the properties of the Gaussian distribution with its 


growth. Despite possible presence of finite-size effects, an 
industry-related clustering structure is observed for both 
exact and approximate couplings. The performed cutoff 
analysis suggests that the biggest positive couplings as 
well as eigenmodes with the biggest eigenvalues contain 
all information about this structure. Scaling properties 
of the couplings between a small random subset of stocks 
suggest that the underlying network structure is also pre¬ 
served with the number of stocks. Finally, the pairwise 
interaction model also allows one for defining the exter¬ 
nal and internal biases which correspond to contribution 
of external fields and couplings respectively. Discrepan¬ 
cies between them might be used as a precursor of immi¬ 
nent financial instabilities and should be explained from 
a financial point of view. These non-trivial market prop¬ 
erties as well as their historical evolution will be studied 
in the future works, where one of the natural extensions 
to incorporate dynamics is to consider the kinetic Ising 
model [in]- 
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TABLE I. List of the companies which stock prices are used for the calculations in the paper. 


Ticker 

Name 

Sector 1 Ticker 

Name 

Sector 

ABT 

Abbott Laboratories 

Hea 

AIG 

American International Group, Inc. 

Fin 

AMGN 

Amgen Inc. 

Hea 

APA 

Apache Corp. 

Bas 

ABC 

Anadarko Petroleum Corp. 

Bas 

AAPL 

Apple Inc. 

Con 

AXP 

American Express Company 

Fin 

BA 

The Boeing Company 

Ind 

BAG 

Bank of America Corp. 

Fin 

BAX 

Baxter International Inc. 

Hea 

BMY 

Bristol-Myers Squibb Company 

Hea 

C 

Citigroup, Inc. 

Fin 

CAT 

Caterpillar Inc. 

Ind 

CELG 

Celgene Corporation 

Hea 

CL 

Colgate-Palmolive Co. 

Con 

CMCSA Comcast Corporation 

Ser 

COP 

ConocoPhillips 

Bas 

COST 

Costco Wholesale Corp. 

Ser 

CSCO 

Cisco Systems, Inc. 

Tec 

CVS 

CVS Caremark Corp. 

Ser 

CVX 

Chevron Corp. 

Bas 

DD 

E. 1. du Pont de Nemours and Co. 

Bas 

DE 

Deere & Company 

Ind 

DELL 

Dell Inc. 

Tec 

DHR 

Danaher Corp. 

Ind 

DIS 

The Walt Disney Company 

Ser 

DOW 

The Dow Chemical Company 

Bas 

EMC 

EMC Corporation 

Tec 

EMR 

Emerson Electric Co. 

Tec 

EOG 

EOG Resources, Inc. 

Bas 

EXC 

Exelon Corp. 

Uti 

F 

Ford Motor Co. 

Con 

GE 

General Electric Company 

Ind 

HAL 

Halliburton Company 

Bas 

HD 

The Home Depot, Inc. 

Ser 

HON 

Honeywell International Inc. 

Ind 

HPQ 

Hewlett-Packard Company 

Tec 

IBM 

International Business Machines Corp. 

Tec 

INTO 

Intel Corp. 

Tec 

JNJ 

Johnson & Johnson 

Hea 

JPM 

JPMorgan Chase & Co. 

Fin 

KO 

The Coca-Cola Company 

Con 

LLY 

Eli Lilly and Company 

Hea 

LOW 

Lowe’s Companies Inc. 

Ser 

MGD 

McDonald’s Corp. 

Ser 

MDT 

Medtronic, Inc. 

Hea 

MMM 

3M Company 

Cng 

MO 

Altria Group Inc. 

Con 

MRK 

Merck & Co. Inc. 

Hea 

MSFT 

Microsoft Gorp. 

Tec 

NKE 

Nike, Inc. 

Con 

ORCL 

Oracle Corporation 

Tec 

OXY 

Occidental Petroleum Corp. 

Bas 

PEP 

PepsiCo, Inc. 

Con 

PEE 

Pfizer Inc. 

Hea 

PG 

The Procter & Gamble Company 

Con 

PNC 

The PNC Financial Services Group Fin 

SLB 

Schlumberger Limited 

Bas 

SO 

Southern Company 

Uti 

T 

AT&T, Inc. 

Tec 

TGT 

Target Corp. 

Ser 

TJX 

The TJX Companies, Inc. 

Ser 

TXN 

Texas Instruments Inc. 

Tec 

UNH 

UnitedHealth Group Incorporated 

Hea 

UNP 

Union Pacific Corp. 

Ser 

USB 

U.S. Bancorp 

Fin 

UTX 

United Technologies Corp. 

Ind 

VZ 

Verizon Gommunications Inc. 

Tec 

WFC 

Wells Fargo & Company 

Fin 

WMT 

Wal-Mart Stores Inc. 

Ser 

XOM 

Exxon Mobil Corp. 

Bas 





Industry sectors are defined as Basic Materials (Bas), Conglomerate (Cng), Consumer Goods (Con), Financial (Fin), 
Healthcare (Hea), Industrial Goods (Ind), Services (Ser), Technology (Tec) and Utilities (Uti). 



